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ABSTRACT 

Context. Angular momentum transport in accretion discs is often believed to be due to magnetohydrodynamic turbulence mediated 
by the magnetorotational instability (MRI). Despite an abundant literature on the MRI, the parameters governing the saturation ampli¬ 
tude of the turbulence are poorly understood and the existence of an asymptotic behavior in the Ohmic diffusion regime is not clearly 
established. 

Aims. We investigate the properties of the turbulent state in the small magnetic Prandtl number limit. Since this is extremely computa¬ 
tionally expensive, we also study the relevance and range of applicability of the most common subgrid scale models for this problem. 
Methods. Unstratified shearing boxes simulations are performed both in the compressible and incompressible limits, with a resolu¬ 
tion up to 800 cells per disc scale height. The latter constitutes the largest resolution ever attained for a simulation of MRI turbulence. 
Different magnetic field geometry and a large range of dimensionless dissipative coefficients are considered. We also systematically 
investigate the interest of using large eddy simulations (LES). 

Results. In the presence of a mean magnetic field threading the domain, angular momentum transport converges to a finite value in 
the small Pm limit. When the mean vertical field amplitude is such that /3, the ratio between the thermal and magnetic pressure, equals 
10 3 , we find a ~ 3.2 x 10 -2 when Pm approaches zero. In the case of a mean toroidal field for which ji = 100, we find a ~ 1.8 x 1CL 2 
in the same limit. Both implicit LES and Chollet-Lesieur closure model reproduces these results for the a parameter and the power 
spectra. A reduction in computational cost of a factor at least 16 (and up to 256) is achieved when using such methods. 

Conclusions. MRI turbulence operates efficiently in the small Pm limit provided there is a mean magnetic field. Implicit LES offers 
a practical and efficient mean of investigation of this regime but should be used with care, particularly in the case of a vertical field. 
Chollet-Lesieur closure model is perfectly suited for simulations done with a spectral code. 
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1. Introduction 

Determining the rate of angular momentum transport in accre¬ 
tion discs is considered as one of the key unsolved astrophysical 
questions. Accretion is believed to be at the origin of the radia¬ 
tion emitted by some of the most luminous sources in the uni¬ 
verse from active galactic nuclei to X-ray binaries and is also a 
major process at work during planet formation in protoplanetary 
discs (Frank et al. 2002). In addition, accretion discs are ubiq¬ 
uitous in the universe and affect the dynamics, evolution and 
appearance of multiple astrophysical objects at all spatial and 
energy scales. The accretion rate can be indirectly constrained 
by the luminosity of high-energy sources or the life-time of pro¬ 
toplanetary discs, and, if the transport process is modeled by an 
a viscosity (Shakura & Sunyaev 1973), such an estimate gives 
0.01 < a < 0.4 depending on the system considered (King et al. 
2007). 

The physical origin of angular momentum transport has to 
be understood to explain such an efficient radial angular mo¬ 
mentum transport. Currently, the most widely accepted mech¬ 
anism is magnetohydrodynamic (MHD) turbulence induced by 
the non-linear evolution of the magneto-rotational instability 
(MRI, Balbus & Hawley 1991). That instability along with its 
nonlinear development has been extensively studied in the last 


two decades (Balbus & Hawley 1998; Balbus 2003; Fromang 
2013). Using local simulations performed in the framework of 
the shearing box, Hawley et al. (1995) quickly established that 
the MRI develops into vigorous MHD turbulence that efficiently 
transports angular momentum outward, a result that was later 
confirmed to be independent of the field geometry (Hawley et al. 
1996) or to the background disc stratification (Brandenburg et al. 
1995; Stone et al. 1996). Only recently the sensitivity of MRI- 
driven MHD turbulence saturation to small scale dissipation in 
such idealized simulations has been identified: when magnetic 
field diffusion is dominated by an ohmic resistivity ?/, a is an 
increasing function of the magnetic Prandtl number Pm, the ra¬ 
tio between the kinematic viscosity v and //. This is the case 
both in the presence of a mean vertical magnetic field (Lesur & 
Longaretti 2007) and of a mean azimuthal magnetic field (Simon 
& Hawley 2009). Such a “Pm-effecf ’ has a significant impact on 
the rate of angular momentum transport measured in homoge¬ 
neous shearing boxes simulations: in the case of a mean vertical 
field such that the plasma [5 parameter (defined in section 2.4) 
amounts to 10 3 , Longaretti & Lesur (2010) showed that a varies 
by a factor of about 5 when Pm only varies between 1 /4 and 4, 
without any sign of the relation flattening at either range. The 
dynamo case (i.e. no mean magnetic field threading the compu¬ 
tational domain) also displays a large sensitivity to Pm (Fromang 
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et al. 2007; Simon & Hawley 2009; Guan et al. 2009; Bodo et al. 
2011) that persists when density stratification is included (Simon 
et al. 2011). In this dynamo regime and without stratification, the 
effect is even stronger and dynamo action is suppressed for Pm 
values smaller than unity, in which case the flow remains lam¬ 
inar, i.e. a goes to zero. A clear understanding of the origin of 
the effect of Pm on the rate of angular momentum transport is 
still lacking and is currently a matter of active research. The first 
results, based on methods borrowed from the fluid community 
(Herault et al. 2011; Riols et al. 2013) are promising (Riols et al. 
2015) and should be extended to more realistic geometries and 
dimensionless numbers. 

Taken together, the results described above question the rele¬ 
vance of MRI-driven MHD turbulence as the dominant transport 
mechanism in accretion discs where the Prandtl number can be 
orders of magnitude lower (Pm « 1) than has been explored in 
published simulations (Balbus & Henri 2008). There is clearly 
the possibility that a becomes vanishingly small as Pm decreases 
to small, but astrophy sic ally relevant values. The first goal of this 
paper is to investigate that asymptotic behavior by means of high 
resolution simulations performed in the homogeneous shearing 
box. In doing so, we will leave aside the dynamo case and focus 
on vertical and azimuthal mean field configurations. Using high 
resolution simulations (such that the coverage in Pm now ex¬ 
tends over more than two orders of magnitudes, from 10 2 to 4), 
we will show that a asymptotically converges to a well defined 
finite value in both cases. However, the computational cost asso¬ 
ciated with such simulations is extremely large (for example, 15 
millions CPU hours are needed to complete the 1000 3 simulation 
we describe in section 3.1.1 on a BlueGen/Q machine ranked 
42nd on the Top500 supercomputing website in november 2014 
1 ). In practice, such a large cost is prohibitive if one wants to 
perform a parameter survey in the asymptotic regime with addi¬ 
tional physics included and simply prevents global simulations 
to be performed in that limit. The second goal of this paper is 
thus to investigate the possibility to use sub-grid scale models as 
a mean to reduce the computational cost associated with MHD 
turbulence simulations in the small Pm limit. 

Recent years have seen significant progress in our under¬ 
standing of the consequences of ambipolar diffusion (Bai & 
Stone 2011, 2013; Simon et al. 2013) and the Hall effect (Kunz 
& Lesur 2013; Lesur et al. 2014; Bai 2015), both of which are 
particularly important in shaping the structure of protoplanetary 
discs. In this paper we will restrict our attention to Ohmic re¬ 
sistivity as the sole source of magnetic field dissipation. Such a 
regime is relevant to describe the very inner parts of protoplan¬ 
etary discs (at stellar distances of a few tens of an AU) but also 
cataclysmic variables (CV) discs and the outer parts of X-ray 
binaries discs (Balbus & Henri 2008). In principle, the analysis 
we present here should also be carried for such cases where am¬ 
bipolar diffusion or the Hall effect are the dominant magnetic 
field diffusion processes. 

The plan of the paper is as follows. In the following sec¬ 
tion, the equations and numerical methods are given. Section 3 
presents our results, focusing first on simulations that explicitly 
resolve the small dissipative scales (section 3.1) and then on two 
different methods to perform Large Eddy Simulations (LES) in 
section 3.2. We finally conclude and discuss the implications of 
our work in section 4. 


1 see http://www.top500.org/list/2014/ll/ 


2. Methods 

2.1. Equations and notations 

In this paper, we use the shearing box approximation (Hawley 
et al. 1995). Namely, we compute the evolution of the fluid in 
a small box centered at a radius tq of the disc and rotating at 
the same velocity as the fluid at r — n.t- As the size of the box 
is small compared to the radial position, the curvature terms of 
the standard fluid equations can be simplified (Hawley & Balbus 
1992). We thus use Cartesian coordinates ( x,y,z ) with units vec¬ 
tors (/,/, k). In this coordinate system, we note (L x ,L y ,L z ) the 
size of the computational box. Neglecting the vertical compo¬ 
nent of the gravitational acceleration, we solve the following set 
of equations: 


d,p + V-(pv) = 0 (1) 

d t (pv) + V • (pvv - BB ) + V P ro , = 2qpQr 0 xi - 2pQo x v 

+V • T (2) 

d,B = V x (v x B - 77V x B) ( 3 ) 

where p is the fluid density, v its velocity in the rotating frame, 
B is the magnetic field and ih, is the angular velocity of the 
fluid at /'o, corresponding to the angular velocity of the box. q 
stands for the background keplerian shear and is taken equal to 
1.5 throughout this paper. P to , is the total pressure, the sum of 
the thermal pressure P and the magnetic pressure B 2 / 2. Explicit 
dissipation is accounted for through the Ohmic resistivity // and 
the kinematic viscosity v that enters in the viscous stress tensor 
T defined as 

T u = pv idjVi + djVj - ■ vj . (4) 


The amplitude of viscosity and resistivity are set using the mag¬ 
nitude of the Reynolds number Re and magnetic Reynolds num¬ 
ber Rm with the following relations 


Re = 



OqL 2 

n 


which can also serve as an alternative definition for the magnetic 
Prandtl number Pm already given in the introduction: 


v 

Pm = — 
>1 


Rm 

~Re' 


(5) 


We performed simulations that consider two different flavors 
of the system of equations (1), (2) and (3). First, we solve the 
above equations in the incompressible limit. In this case, the den¬ 
sity is constant and equation (1) reduces to V*v = 0. We used the 
code SNOOPY in that case (see section 2.2). In the second type 
of simulations, we solve the full set of equations (and refer to that 
case as the compressible simulations) using the code RAMSES 
(section 2.3). We briefly describe below the two codes and the 
specificities of each set of simulations. 


2.2. Incompressible simulations: the SNOOPY code 

The SNOOPY code is a pseudo spectral code that solves the 
incompressible MHD equations in a Fourier basis. It uses a low- 
storage 3 rd order Runge-Kutta integrator and works in a sheared 
frame comoving with the mean flow, which is equivalent to a 
3 rd order Fargo scheme (Masset 2000). SNOOPY uses a 3/2 
antialiasing rule to eliminate the excitation of spurious modes 
during the computation of quadratic nonlinearities. 
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The use of a sheared Fourier basis makes the boundary con¬ 
ditions periodic in the y and z directions and shear-periodic in 
the a - direction. SNOOPY conserves linear momentum and mag¬ 
netic flux down to machine precision. Since spectral methods are 
inherently diffusion-free and energy conserving, the addition of 
dissipation is required to mimic the damping due to small scale 
dissipation processes 2 . In SNOOPY, one can choose second or¬ 
der diffusion operators, hyperdiffusion operators or Chollet & 
Lesieur (1981) subgrid models (see section 3.2.1). 

2.3. Compressible simulations: the RAMSES code 

The RAMSES code is a finite volume code that solves the com¬ 
pressible MHD equations on a Cartesian grid (Teyssier 2002; 
Fromang et al. 2006) using the constrained transport algorithm 
(Evans & Hawley 1988). We use a version of the code for which 
the grid is uniform (i.e. without the adaptive mesh refinement). 
The source terms associated with the tidal potential are included 
as described by Stone & Gardiner (2010). We use the HLLD 
Riemann solver (Miyoshi & Kusano 2005) with monotonized 
central slope limiter, shearing box boundary conditions in the x 
direction (Hawley et al. 1995) and periodic boundary conditions 
in the y and z directions. For extended box sizes such as used in 
the mean vertical field case (see section 2.4), it has been found 
that radially variable numerical dissipation causes the turbulent 
stress to vary accross the box. To avoid that problem, we used 
the FARGO algorithm (Masset 2000; Stone & Gardiner 2010) in 
that case, which also improves the efficiency of the code through 
an increase in the timestep. 

Throughout this paper, we use an isothermal equation of 
state to close the system of MHD equations in compressible sim¬ 
ulations, such that P = pcQ, where co stands for the constant 
sound speed. We start the simulation with an initial uniform den¬ 
sity p = po and choose L z — H = cq/DIq where H is the disc scale 
height. 

2.4. Models parameters 

We start the simulations with a uniform initial magnetic field: 

B = Boyj + Bo z k. ( 6 ) 

Two different initial magnetic configurations are considered in 
the following, namely pure azimuthal field (B t)z = 0) simulations 
and pure vertical field simulations (Boy = 0). The strength of the 
magnetic field is defined using the plasma parameter /?, that is 
given by 


( rpoc-Q 

«. = J 4 , 

Pi ~ ) (qClL.J 2 

l 4 

in compressible runs, 

in incompressible runs. 

(7) 



(8) 


For the mean azimuthal (resp. vertical) magnetic field simula¬ 
tions, we used /3 y = 112.5 (resp. /3 Z = 10 3 ). At the start of each 
runs, small amplitude random velocity perturbations are added 
on the three velocity components with an amplitude equal to 1 % 

2 Note that the numerical stability of the scheme does not require 
any dissipation. However, the absence of dissipation in a turbulent flow 
naturally leads to thermalisation, a situation which does not occur in 
natural systems which always exhibit dissipation processes. Numerical 
dissipation is therefore required on physical grounds to break the ther¬ 
modynamic equilibrium and create the well known energy cascade pic¬ 
ture. 


of the sound speed in compressible runs. In incompressible sim¬ 
ulations both velocity and magnetic random perturbations are 
added with an amplitude of 0.1 QL-. 

The size of the computational box is either (L Z ,4L Z , L z ) 
for the simulations with a mean azimuthal magnetic field, or 
(4L Z , 4L Z , L z ) for the simulations with a mean vertical magnetic 
field. In the latter case, it is indeed well known that boxes with 
L x = L- artificially enhance the importance of recurrent bursts 
in the flow structure (Bodo et al. 2008; lohansen et al. 2009). 
In such box sizes and with such vertical magnetic field, Bai & 
Stone (2014) recently reported zonal flows. We also found such 
structures in our simulations. 

The resolution varies from 32 cells per unit length up to 800 
cells per unit length. Such large resolutions allow us to reach the 
largest Reynolds number (run Y-C-Re85000, Re = 85000) and 
the smallest Prandtl number (run Z-I-Re40000, Pm — 0.01) ever 
published. As larger structures are expected in the y-dircction, 
we typically use a resolution that is twice as coarse in that direc¬ 
tion. 


3. Results 


The whole set of runs discussed in the remaining of this paper 
is listed in tables 1, 2 and 3, where the first column provides 
the simulation labels. Runs starting with a pure azimuthal (verti¬ 
cal) magnetic field are labelled with letter “Y” (“Z”). Likewise, 
compressible (incompressible) simulations are labelled with let¬ 
ter “C” (“I”). Large Eddy Simulations (LES) are labelled either 
“ILES” or “CL” depending on the subgrid scale model that is 
used (see section 3.2). The remaining columns in tables 1, 2 and 
3 give the run resolution (col. 2) and duration Tr u „ (col. 3), the 
Reynolds number (col. 4), the magnetic Reynolds number (col. 
5), the magnetic Prandtl number (col. 6), and a (col. 7), the sum 
of aR ey (col. 8) and ctMax (col. 9). The later are defined by the 
following relations: 


OtRey 


& Max 


{pSv r Sup 
Po 

<SVrSup 
(ilLz ) 2 
(B r Bp 

<£,%> 

Po(HLz ) 2 


in compressible runs, 
in incompressible runs, 

in compressible runs, 
in incompressible runs. 


(9) 

( 10 ) 


Here (.) denotes an average over the shearing box volume and 
over time and <Sv is the velocity difference to the laminar sheared 
flow. Except for the shortest runs (see for instance section 3.1.1), 
the turbulent transport rates as measured by the a parameters are 
time averaged over the last 60 orbits of the models. The last col¬ 
umn finally gives the ratio between magnetic and hydrodynamic 
transport rate. 


3.1. Direct Numerical Simulations (DNS) in the small Pm 
regime 

In this section, we present the results of the resolved simula¬ 
tions (meaning that kinematic viscosity and ohmic resistivity are 
explicitly included in the calculation), focusing on the rate of 
angular momentum transport and on the power spectra of the 
turbulent flow. 


3.1.1. A 1000 3 MRI simulation 

Because of the tremendous computational cost that was associ¬ 
ated with that simulation, we start with a description of model 
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Model 

Resolution 

T R, m 

Re 

Rm 

Pm 

a 

&Mcix 

&Rey 

R 

Y-C-Re650 

(64,128,64) 

100 

650 

2600 

4 

3.0 x 10^ 2 

2.4 x 10~ 2 

5.9 x 10~ 3 

4.2 

Y-C-Re2600 

(64,128,64) 

100 

2600 

2600 

1 

2.5 x 10~ 2 

1.9 x 10~ 2 

5.7 x 10- 3 

3.4 

Y-C-Re 13000 

(128,256, 128) 

100 

13000 

2600 

0.2 

1.8 x 10~ 2 

1.4 x 10~ 2 

4.6 x 10 3 

3.0 

Y-C-Re26000 

(256,512,256) 

100 

26000 

2600 

0.1 

2.0 x 10^ 2 

1.5 x 10~ 2 

5.0 x 10’ 3 

3.0 

Y-C-Re85000 

(800,1600,832) 

35 

85000 

2600 

0.03 

1.8 x 10~ 2 

1.3 x 10~ 2 

4.6 x 10- 3 

2.9 

Y-ILES-C-64 

(64,128,64) 

100 

- 

2600 

- 

1.9 x 10~ 2 

1.4 x 10~ 2 

4.7 x 10’ 2 

3.0 

Y-ILES-C-128 

(128,256, 128) 

100 

- 

2600 

- 

1.9 x 10- 2 

1.4 x 10~ 2 

4.8 x 10 3 

3.0 

Y-ILES-C-256 

(256,512,256) 

100 

- 

2600 

- 

1.9 x 10~ 2 

1.4 x 10~ 2 

5.0 x 10- 3 

2.8 

1SES runs with a mean azimuthal field. The box size is (L x , L y 

,4) = 

(4,4L-, 4). 

T nm is given in local orbits time ui 


gives for each run its resolution, duration, Reynolds number, magnetic Reynolds number, magnetic Prandtl number, total stress. Maxwell stress, 
Reynolds stress and the ratio of the two. 


Model 

Resolution 

T Run 

Re 

Rm 

Pm 

a 

&Max 

&Rey 

R 

Z-C-Re400 

(128,64,32) 

100 

400 

400 

1 

6.5 x 10~ 2 

4.0 x 10~ 2 

2.4 x 10^ 2 

1.6 

Z-C-Re800 

(256,128,64) 

100 

800 

400 

0.5 

5.1 x 10~ 2 

2.9 x 10~ 2 

2.2 x 10~ 2 

1.3 

Z-C-Re3000 

(512,256,128) 

100 

3000 

400 

0.13 

3.7 x 10~ 2 

2.0 x 10~ 2 

1.8 x 10^ 2 

1.1 

Z-C-Re8000 

(1024.512,256) 

100 

8000 

400 

0.05 

3.1 x 10~ 2 

1.5 x 10~ 2 

1.6 x 10^ 2 

1.0 

Z-ILES-C-32 

(128,64,32) 

100 

- 

400 

- 

4.1 x 10~ 2 

2.2 x 10- 2 

1.9 x 10~ 2 

1.2 

Z-ILES-C-64 

(256,128,64) 

100 

- 

400 

- 

3.6 x 10~ 2 

1.8 x 10~ 2 

1.8 x 10~ 2 

1.0 

Z-ILES-C-128 

(512,256,128) 

100 

- 

400 

- 

3.3 x 10~ 2 

1.6 x 10~ 2 

1.7 x 10^ 2 

1.0 

Z-ILES-C-256* 

(1024,512,256) 

40 

- 

400 

- 

3.3 x 10~ 2 

1.6 x 10- 2 

1.8 x 10~ 2 

0.9 


Table 2. Same than table 1 but for RAMSES runs with a mean vertical field. The box size is (L x , L v , L z ) = (4L-, 4L-, L z ). *To save computational 
time, this model was performed by restarting model Z-C-Re8000 at t = 50 removing explicit viscosity 


Model 

Resolution 

T«„„ 

Re 

Rm 

Pm 

a 

& Max 

&Rey 

R 

Z-I-Rel300 

(256,256,64) 

53 

1333 

400 

0.3 

3.8 x 10~ 2 

2.8 x 10 2 

9.9 x 10~ 3 

2.8 

Z-I-Re20000 

(1536,768,192) 

53 

20000 

400 

0.02 

3.5 x 10~ 2 

2.6 x 10 2 

9.2 x 10~ 3 

2.8 

Z-I-Re40000 

(3072,1536,384) 

53 

40000 

400 

0.01 

3.3 x 10’ 2 

2.3 x 10’ 2 

1.0 x 10~ 2 

2.4 

Z-CL-192-a 

(768,384,192) 

53 

- 

400 

- 

3.3 x 10~ 2 

2.4 x 10- 2 

8.8 x 10~ 3 

2.8 

Z-CL-192-b 

(768,384,192) 

53 

- 

400 

- 

3.2 x 10~ 2 

2.3 x 10~ 2 

9.3 x 10~ 3 

2.5 

Z-CL-128-a 

(512,256,128) 

53 

- 

400 

- 

3.5 x 10’ 2 

2.4 x 10’ 2 

1.0 x 10~ 2 

2.4 

Z-CL-128-b 

(512,256,128) 

53 

- 

400 

- 

3.4 x 10’ 2 

2.5 x 10 2 

9.7 x 10~ 3 

2.5 

Z-CL-96-a 

(384,192,96) 

53 

- 

400 

- 

3.4 x 10~ 2 

2.5 x 10~ 2 

9.2 x 10~ 3 

2.8 

Z-CL-96-b 

(384,192,96) 

53 

- 

400 

- 

3.7 x 10’ 2 

2.6 x 10’ 2 

1.1 x 10~ 2 

2.3 


Table 3. Same than table 1 but for SNOOPY runs. All runs have a box size of ( L x , L,,, L z ) - (4L z , 4L z , L-). 



Y-C-Re85000, performed with RAMSES. In this model. Re = 
85000 and Pm = 0.03. This simulation was performed with a 
resolution (N x , N y , N z ) — (800,1600,800). Based on past ex¬ 
perience and published results of simulations using the same 
setup at lower Reynolds number with a similar code (Simon 


& Hawley 2009), we are confident that the small scales of the 
flow are properly accounted for in that simulation. In the x and 
z directions, the number of cells thus amounts to 800 cells per 
disc scaleheight, which is the highest resolution ever achieved 
of MRI-driven MHD turbulence with a 2nd order compressible 
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Fig. 1. Top panel: Time evolution of a in model Y-C-Re85000 for which 
Pm = 0.03. The initial 15 orbits were performed in the ideal MHD 
regime, before including explicit diffusion coefficients. The hatched re¬ 
gion represents the time interval over which the angular momentum 
transport is averaged. Bottom panel: Y-C-Re85000 run, kinetic (solid 
green line) and magnetic energy (dashed blue line) power spectrum. 


code. This gigantic simulation was run at the IDRIS supercom¬ 
puting center on the BlueGen-Q machine Turing, using 32768 
cores. We found that the most efficient configuration was to use 
4 threads per cores, meaning that a total of 131 072 threads (or, 
equivalently, MPI sub-domains) were used. We used approxi¬ 
mately 10 hours of CPU time per timesteps. Running the sim¬ 
ulation for over 35 orbits (~ 10 6 timesteps) thus required about 
10 7 hours of CPU time, corresponding to 300 hours of wall clock 
time (i.e. about two weeks). 

As described by Simon & Hawley (2009), in the presence 
of a mean azimuthal magnetic field, resistivity can prevent the 
linear instability from transiting to a turbulent state when the 
simulation starts from a laminar flow, even though it is found 
to remain turbulent on long timescales when that linear phase 
is bypassed. We thus adopted the method described by these au¬ 
thors and performed the run is two steps: we first solved the ideal 
MHD equations (i.e. with vanishing viscosity and resistivity pa¬ 
rameters). The only dissipative effects are numerical in origin 
during that part of the calculation. A turbulent state is reached 
after ~ 10 orbits. The turbulent transport as measured by a dis¬ 
plays fluctuations around a well defined mean value of about 
6x10 2 for the next few orbits (see figure 1, top panel). At 1 - 15, 
we restarted the simulation for an additional 20 orbits with dissi¬ 
pation coefficients such that Re = 85000 and Pm = 0.03. As seen 
on the top panel of figure 1, the turbulence amplitude is quickly 
modified by the presence of those dissipation coefficients and 


reaches a new steady state after a few additional orbits 3 . The 
nature of the flow is illustrated by a series of snapshots of the 
last time step of the run in figure 2. On the density plot one can 
recognize large scales density waves and low amplitude shocks. 
The turbulent magnetic field is dominated by large scale struc¬ 
tures. On the contrary, the velocity snapshot shows both large 
and small scales structures, as expected given the very small Pm 
used here. 

We next averaged the transport rate after the system has 
reached a quasi steady state and until the end of the simulation. 
This period corresponds to t > 21 orbits and is hatched on 
figure 1 (top panel). The value of aR ey , ctMax and a we obtained 
are 4.6 x 1(T 3 , 1.3 x 1CU 2 and 1.8 x 10 -2 , respectively (see 
also Table 1). This a value is comparable to the one reported 
by Simon & Hawley (2009) for their most resolved run, for 
which Pm = 0.25, Rm = 3200 and (i - 450. The Maxwell 
and Reynolds stresses display a ratio of ~ 3 that is typical of 
such simulations. There is of course a degree of arbitrariness 
regarding the exact time at which we decide that the system 
has reached “a quasi steady state” and start averaging. We have 
checked that the statistics we consider do not vary significantly 
when making modest changes to the averaging period. For 
example, we find a = 1.9 x 10~ 2 and a = 1.8 x 10~ 2 when 
averaging over the periods t > 18 and t > 25 orbits, respectively. 
This is only a variation by 7% and give an idea of the uncertainty 
associated with that measurement. 

To obtain a more accurate view of the energy budget at each 
scales, we also show on figure 1 (bottom panel) the time av¬ 
eraged power spectra. Because of the shearing box boundary 
conditions, we consider time dependent unsheared wave vec¬ 
tor k and the shell filter decomposition of the physical variables 
(Hawley et al. 1995). Being spherically symmetric, we note that 
this method filters out the information about the flow anisotropy 
which is known to be significant for MRI-driven turbulence (see 
Murphy & Pessah 2015, Figure 3. of Lesur & Longaretti 2011 
and section 3.2.1). For each wavenumbers, the kinetic and mag¬ 
netic energy are then given by 


E K (k) = v 2 (k)/2 (11) 

E„(k) = B 2 k (k)/ 2, (12) 

respectively, where the bar denotes an average over time. 
Magnetic energy dominates at large scales ( k' < 15, where 
we have defined k' — k/ln so that it corresponds to scales 
l' > H/ 15) while kinetic energy dominates at smaller scales. 
This is a signature of the small Pm value of the simulation, which 
implies that the resistive dissipation length is larger than the 
viscous dissipation length (Schekochihin et al. 2004). At small 
scales, the motions essentially correspond to hydrodynamic tur¬ 
bulence associated with a forward cascade (Lesur & Longaretti 
2011). We will exploit that scale separation in section 3.2 when 
designing sub-grid scale models for the hydrodynamic part of 
the flow. The kinetic energy power spectrum follows a power 
law with exponent -3/2 over the range 2 < k! < 20, thus cover¬ 
ing one order in magnitude in spatial scales. The same expo¬ 
nent has been reported recently in other high resolution sim¬ 
ulations of MRI-driven MHD turbulence both in the dynamo 
regime (Fromang 2010) and in the presence of a net vertical 

3 Note that all the runs with a mean azimuthal magnetic field are 
executed using the same procedure. However, except for model Y-C- 
Re85000 for which the computational cost is considerable, they have 
usually been run for 100 orbits to obtain a more precise measure of a. 
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Fig. 3. Mean value of the angular momentum transport (measured by means of a) for the mean toroidal field simulations (blue circles ), the mean 
vertical field simulations performed with RAMSES (green filled circles) and with SNOOPY (green filled squares). The dashed lines are power-law 
functions that approximately describe the data (see text). 


field (Lesur & Longaretti 2011), suggesting it is a general feature 
of the flow. We note that homogeneous forced MHD turbulence 
also displays the same exponent (Mason et al. 2008), although 
the result is still debated (Beresnyak 2014). The surprising re¬ 
sult here is that the magnetic energy does not show any obvious 
signature of a power-law regime, contrary to the case of homo¬ 
geneous and driven MHD turbulence, while still being compa¬ 
rable in magnitude with Eg. More work is needed to clarify the 
reason for that discrepancy and to better understand the origin 
of the -3/2 exponent that is obtained in the case of MRI-driven 
MHD turbulence. 

3.1.2. Angular momentum transport 

We plot on figure 3 the total averaged stress a for all the resolved 
simulations we performed. Error bars cr a are estimated with the 
method presented in Longaretti & Lesur (2010). In general, we 
find cr a ~ 5 x 10 for the vertical field model (thus giving 
(Tala ~ 5%). This is consistent with the estimate of Longaretti & 
Lesur (2010). The error estimate is smaller in the azimuthal field 
case, for which we obtained cr a ~ 10 3 which corresponds to 
(Ta/a ~ 1%. We note that we could not apply the same method 
for model ’Y-C-Re85000’ because of the short duration of the 
integration in that case. A simple standard deviation is thus plot¬ 
ted instead and explains the larger error bar for that particular 
model. 

In agreement with previous results (see section 1) we find 
that a increases with Pm. However, the main result of figure 3 


(and the main result of the paper) is that there is now convincing 
evidence of the convergence of the angular momentum transport 
rate at small Pm toward a well defined, nonzero value. This is 
the case for both field geometry. 

The azimuthal field case: In this case, and with constant Rm 
and p, we find that a fit to the data is given by the formula 

a By = dByrnin + CyPftl °' 58 , (13) 

with the values 4 of asymin = 1-8 x ICR 2 and C y = 5.5 x ICE 3 , 
respectively. This means that a good estimate of a is already 
obtained at Pm = 0.2 for which a resolution of 128 cells per unit 
length is sufficient. Indeed, we obtained a — 1.8 x ICE 2 which is 
equal to the asymptotic value of a at vanishingly low Pm. 

The vertical field case: Here, we find a somewhat larger trans¬ 
port coefficient that can be fitted by the relation 

&Bz — tT Bzniin T C-Plll (14) 

with a b zmin = 3.2 x ICE 2 and C, = 3.3 X ICE 2 , with fixed Rm and 
p. Moreover, there is good agreement between the compressible 
and incompressible approach in this vertical field case: for both 
flow type the simulations converge at small Pm toward the same 

4 Note that the asymptotic transport values a Bym i„ and a Bzmin depend 
in principle on Rm and p. See the discussion in section 4. 
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Fig. 4. Left panel: Time history of a (blue curve), a Max (red cur\’e) and a Rey (green cur\’e) in model Z-C-Re3000. Right panel: scatter plot showing 
the total angular momentum transport rate a as a function of a axi which measures the angular momentum transport due to k' = 0 and k' z = 1 modes 
for 120 dumps evenly spaced over model Z-C-Re3000. The green curve plots an approximate fit to the data (see text for details). 


a. As for the azimuthal field case, a good estimate of the asymp¬ 
totic rate of angular momentum transport is already obtained for 
Pm = 0.13 (using a resolution of 128 cells per unit length), for 
which we found a = 3.7 x 10 2 , i.e. a value that differs by about 
15% from the asymptotic limit. 

As already noted (section 2.4), in the presence of a vertical 
magnetic field we find that the time history of a displays sig¬ 
nificant variability. This is illustrated on figure 4 (left panel) for 
the particular case of model Z-C-Re3000 (Pm - 0.13): while 
the time averaged transport rate amounts to a - 3.7 x 10 2 in 
that case, there are numerous peaks during which it reaches val¬ 
ues as high as 0.1 that occur with a typical period of about 5 to 
10 orbits. Such bursts are not unheard of in unstratified shear¬ 
ing boxes with a mean vertical field (Bodo et al. 2008; Latter 
et al. 2009) and have been attributed to recurrent ‘channel-like’ 
modes associated with the MRI. For the set of parameters we 
have considered here, the time history of a suggests that they 
contribute significantly to the turbulent transport. In an attempt 
to quantify that contribution, we have calculated for model Z- 
C-Re3000 the value a axl of the transport that is due to axisym- 
metric ‘channel-like’ modes (for which k' y - 0 and k' - 1) for 
120 dumps evenly spaced between t - 40 and t = 100. The 
mean value of a ax j, averaged over all the dumps of the simula¬ 
tion, amounts to 1.2 x 10 2 : this means that axisymmetric modes 
with k( = 1 accounts for ~ 30% of the angular momentum trans¬ 
port. In agreement with the results of Longaretti & Lesur (2010), 
angular momentum transport is dominated on average by non- 
axisymmetric modes even if the contribution of ‘channel-like’ 
modes is significant. The scatter plot showing the relation be¬ 
tween a and a axi for those 120 dumps is shown on the right 
panel of figure 4, along with an indicative fit of the data given 
by 


a — a o 


/ \ 0 - 6 
/ a axi \ 

\0!axio) 


(15) 


with ofo = 10 2 and a ax m = 10 5 . The positive correlation be¬ 
tween a and a axl indicates that the relative contribution of the 
transport mediated by axisymmetric modes with k( - 1 increases 
during the bursts of activity: for example, a ax i amounts to only 
about 10% of the transport when a - 10 2 but, as indicated by 



Fig. 5. Time history of a Max (blue cur\’e) and a Rey (green curve) in 
model Z-C-Re3000. 


Eq. (15), can contribute for up to 50% of the turbulent activity 
when a - 0.1. These variations are consistent with the results 
of Latter et al. (2009) and with the idea that the bursts seen on 
figure 4 are due to large scale ‘channel-like’ modes (Bodo et al. 
2008). We have repeated the same analysis for all the models 
and we have found a weak dependance of the relative fraction 
of axisymmetric transport with Pm: a ax i/a = 0.39, 0.34, 0.32 
and 0.29, respectively for Pm = 1, 0.5, 0.13 and 0.05 and simi¬ 
lar results for the incompressible runs with a ax i/a - 0.38, 0.35, 
respectively for Pm - 0.02 and 0.01, 

It is also noteworthy that the angular momentum transport 
in the presence of a vertical field is equally due to Maxwell and 
Reynolds stresses in the compressible simulation, whereas the 
classical ratio of approximately 3 is obtained in an incompress¬ 
ible fluid or with an azimuthal field configuration for the specific 
values of /? and Rm chosen in our investigation. The exact value 
of Maxwell to Reynolds stress ratio (R) are given in the last col¬ 
umn of Tables 1,2, and 3. In an attempt to have a better under¬ 
standing of the relative contribution of Maxwell and Reynolds 
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Fig. 6. Mean value of the angular momentum transport (measured 
by means of a) for the mean toroidal field LES simulation (implicit 
LES blue circles ), the mean vertical field LES simulation performed 
with RAMSES (implicit LES, green filled circles) and with SNOOPY 
(Chollet-Lesieur LES, green filled squares). The resolution is given in 
number of cells per scaleheight. The dashed lines represents the asymp¬ 
totic limits obtained with DNS. 


stresses, we plot on figure 5 their time history over two burst 
cycles for the particular case of model Z-C-Re3000. This plot 
reveals that the Reynolds stress is larger than the Maxwell stress 
during the decaying part of the burst when the ‘channel-like’ 
modes (with k y = 0 and k z = 1) are destroyed by the non-linear 
turbulent dynamics. The computation of the contribution of the 
axisymmetric modes shows that the Maxwell transport is dom¬ 
inated by ‘channel-like’ modes whereas the Reynolds transport 
is due to non-axisymmetric modes. Because of the difference 
between the incompressible and compressible runs, we further 
speculate that such a destruction is associated with the excitation 
of compressible modes such as density waves. However this de¬ 
tailed study, not directly related to the angular momentum trans¬ 
port rate, is beyond the scope of this paper. 

3.2. Large Eddy Simulations in the small Pm limit 

The previous results have shown that angular momentum trans¬ 
port converges toward a well defined limit at small P m . However, 
such simulations are very computationally expensive. Here, we 
investigate the possibility to use a sub-grid scale model instead 
of standard kinematic viscosity. The aim is to reduce the com¬ 
putational cost of the simulations without compromising the 
physics one may want to consider, such as the accretion rate or 
the turbulent power spectra. 

Historically, the study of small Pm flows has mostly re¬ 
lied on simulations using hyperviscosity such as geodynamo 
models (Glatzmaiers & Roberts 1995) or small scale dynamo 
theory (Schekochihin et al. 2007). However, hyperviscosity is 
known to produce numerous artefacts in hydrodynamic turbu¬ 
lence, such as spectral bottlenecks, reduced intermittency and 
spurious isotropisation (Frisch et al. 2008). For this reason, we 
instead focus on Large Eddy Simulations (LES) in which the 
dissipation adapts dynamically to the turbulent cascade initiated 
by the large scales, thereby limiting the artefacts commonly due 
to hyperviscosity. 


LES are routinely used in industrial applications to model 
flows at high Re. However, their MHD counterpart are not 
widespread, the reason being the higher complexity of the 
MHD turbulent cascade compared to the purely hydrodynamic 
cascade. Several subgrid-MHD models have been discussed 
in the literature including Smagorinsky (1963) type models 
(Yoshizawa 1987) which can be used in finite difference/volume 
codes and Chollet & Lesieur (1981) type models (Baerenzung 
et al. 2008) targeted to spectral methods. However, most of these 
models are still under development and their applicability to 
Pm <K 1 flows is yet to be proven. 

In this work we have chosen to use well-known hydrody- 
namical subgrid models to treat the kinetic turbulent cascade 
only, leaving the induction equation with standard Ohmic resis¬ 
tivity. This approach is valid provided that the subgrid model is 
introduced at a scale much smaller than the resistive scale, so 
that the energy cascade is mostly hydrodynamical. Moreover, it 
has the advantage of using well-tested subgrid models which are 
reasonably simple to implement numerically. This type of meth¬ 
ods has already been used to study Taylor-Green flows (Ponty 
et al. 2004) and dynamo action (Ponty et al. 2005) in the limit 
Pm —> 0. We therefore reproduce this approach using very simi¬ 
lar tools in the MRI turbulence context. 

3.2.1. Chollet-Lesieur model in incompressible simulations 

Since our incompressible simulations are done with a spectral 
code, we have used the spectral subgrid model of Chollet & 
Lesieur (1981) (CL) to perform our incompressible LES. This 
model replaces the standard viscosity v by the following expres¬ 
sion in Fourier space: 

vcm = m (0.267 + 9.21 exp[-3.03 k c /k]) (16) 

where k c is the cutoff scale and E(k c ) is the kinetic energy at the 
cutoff scale 5 . This expression encloses long range nonlinear in¬ 
teractions via a constant viscosity at k <zr k c and a cusp due to lo¬ 
cal energy transfers close to the cutoff scale k c . Interestingly, k, is 
the only free parameter of this model. The amount of viscosity it¬ 
self is automatically adjusted according to the amount of energy 
at the cutoff scale. It should be emphasized that this expression is 
only valid for 3D homogeneous and isotropic Kolmogorov tur¬ 
bulence. In principle, it is therefore not applicable to turbulent 
shear flows found in shearing box models since such flows are 
anisotropic (see section 3.1.1). As a first attempt to quantify that 
anisotropy, we proceed as follows. We define a decomposition 
in spherical harmonics 6 such that 

E K (k) = \\v(k)\ 2 = Yj c'J(\k\)Y jm (k), (17) 

jm 

and similarly for the magnetic energy. We focus on the / = 2 
coefficients of that decomposition and define: 

m—2 

*2 = Y l C 2 | 2 /|C 0 | 2 - (18) 

m=-2 

5 Note that several forms for the CL viscosity may be found in the 
literature since this expression is a fit to numerical EDQNM (Eddy- 
Damped Quasi-Normal Markovian) calculations. Our expression comes 
from the asymptotic viscosity of Chollet & Lesieur (1981) and a fit close 
to k c of Sagaut (2006). 

6 This use of spherical harmonics to estimate the anisotropy of tur¬ 
bulence is a standard procedure to study sheared flows (e.g. Biferale & 
Vergassola 2001). 
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Fig. 7. Estimate of the anisotropy at each scale. For each wavenumber, 
the energy is decomposed into spherical harmonics; the anisotropy is 
estimated from the sum of the second order coefficients normalized by 
the isotropic coefficient. 


ao takes significant values when there are strong variations of 
the energy over the shell, or, in other words, when the flow dis¬ 
plays anisotropy at that scale. As shown in figure 7, it decreases 
with k but is still high at small scales for the kinetic energy. 
Nevertheless we will see that the Chollet-Lesieur model gives 
satisfactory results. 

We have performed several simulations using CL subgrid 
model, varying k c and the resolution (runs Z-CL-XXX-X). 
Models labelled a have k c /k max = 0.75 and models labelled 
b have k c /k mm = 1 where k max is the maximum accessible 
wavenumber of the simulation (Table 3). As plotted on figure 
6 (green squares), all of our models recover the statistical results 
of DNS with a reduced resolution. These encouraging results 
are confirmed by turbulent spectra compared to DNS simulation 
(figure 8). We find that both kinetic and magnetic spectra agree 
to less than 10% down to k ~ 20 with the highest resolution DNS 
Z-I-Re40000. We also find that simulations with k/k c = 0.75 ex¬ 
hibit the same convergence properties, at least for the resolution 
we studied. Therefore, the cutoff scale does not seem to have 
much impact on these results, provided it is below the resistive 
dissipation scale. These results demonstrate that CL models can 
be used efficiently to study low Pm flows, with a gain in resolu¬ 
tion of at least a factor 2 associated with a gain in computation 
time of at least a factor 20. 

3.2.2. Implicit LES in compressible simulations 

In the small Pm limit, the simplest possible sub-grid scale 
model when using finite volume codes is certainly the so-called 
“Implicit LES” (ILES). The idea is to capture Ohmic resistivity 
explicitly in the simulation (since it occurs at large scale) but let 
numerical dissipation handle kinetic energy dissipation at small 
scales. This method has already been used in some previous 
works (Fleming et al. 2000; Inutsuka & Sano 2005; Okuzumi 
& Hirose 2011; Flock et al. 2012, 2015), but its range of valid¬ 
ity has never been systematically and quantitatively investigated. 
In this section, we compare the results of such simulations, per¬ 
formed with RAMSES at various resolution with the results pre¬ 
sented in section 3.1. 




Fig. 8. Ratio of power spectra of kinetic energy (top) and magnetic 
energy (bottom) obtained with a Chollet-Lesieur subgrid model with 
kjkmux = 1 and full DNS calculation at Re=40000. 


The azimuthal field case: the magnetic Reynolds number is 
here fixed to Rm = 2600 as for the DNS. We performed a se¬ 
ries of simulations with resolutions ranging from 64 to 256 cells 
per scale height. For all cases, we obtain a — 1.8 x 1CL 2 , in very 
good agreement with our best resolved simulations at Pm = 0.03 
(see figure 6). We then compared the kinetic and magnetic en¬ 
ergy power spectra in these simulations with the spectra obtained 
in the DNS with the lowest Prandlt number (Y-C-Re85000). In 
the upper panel (resp. lower panel) of figure 9, we plot the ratio 
of the kinetic (resp. magnetic) power spectra between the ILES 
and the DNS. The runs give an acceptable agreement with the 
DNS at all available scales (except at small scales in the kinetic 
energy, which is a result of the different hydrodynamical dissipa¬ 
tion): at all resolutions, the deviations to the DNS run for k' < 10 
for both the kinetic and the magnetic energy spectra, are at most 
of the order of 15%. 

The vertical field case: the resolution in this case was varied 
from 32 cells per scale height to 256 cells per scale height. As 
can be seen on figure 6, a decreases as resolution increases, most 
likely as a result of the decrease of the effective numerical vis¬ 
cosity, and converges to the value determined by the DNS sim¬ 
ulations. Quite surprisingly though, the convergence as resolu¬ 
tion is increased toward the asymptotic value of a at low Pm 


9 



























Meheut et al.: MRI turbulence in the small Pm limit 



51 

S' 

6 ? 



k/2 ti 



Fig. 9. Ratio between the power spectra obtained with RAMSES in the 
ILES (green, red, blue and dotted cunes) and model Y-C-Re85000 
in the azimuthal magnetic field cases. Top panel is for kinetic energy 
power spectrum and bottom panel for magnetic energy power spectrum. 
On both panels, the black horizontal line mark the location of unity ra¬ 
tio. 


is slower than the azimuthal field case, despite Rm being much 
larger in that case. Indeed, with 64 cells per scale height, the dif¬ 
ference between the two a values is still of order 20% while it 
has reached convergence in the case of a mean B y . This is prob¬ 
ably due to the stress dependence with Pm being steeper in the 
presence of a mean vertical field. In any case, 128 cells per scale 
height are needed to reach the asymptotic a value to less than 
10%. We note that this is still a gain of a factor of 16 in com¬ 
puting time. The power spectra displayed in figure 10 confirm 
that result: there is a difference of about 50% between the power 
spectra of the ILES and the DNS at large scale for a resolution 
of 32 cells per scale height. Even the largest resolution reveals 
differences of up to 20% in both the kinetic and magnetic power 
spectra, even if the largest scales are converged to better than 
10%, as expected from the good agreement between the trans¬ 
port coefficients in both simulations (figure 6). 

4. Conclusions 

We briefly summarise the main results and ideas of the paper and 
discuss some of their limitations. 

Our main goal was to investigate the properties of the turbu¬ 
lent state in the small magnetic Prandtl number limit by mean 



k/2it 


Fig. 10. Same as figure 9, but for the vertical field case. 


of local DNS simulations. We showed that in the presence of a 
mean magnetic field threading the domain, angular momentum 
transport converges to a finite value in the small Pm limit. This 
result is valid both with a vertical and with an azimuthal mean 
magnetic field with the asymptotic values at small Pm being re¬ 
spectively: a — 1.8 X 10 2 with Rm = 2600 and J3 y ~ 10 2 and 
a — 3.2 x 10 2 with Rm = 400 and J3 Z = 10 3 . The obtained val¬ 
ues with our set of parameters are consistent with the estimations 
computed from the lifetime and accretion rate of protoplanetary 
discs (a is a few 10~ 2 ). Obviously, a word of care is in order 
here: the magnetization of accretion discs, such as protoplane¬ 
tary discs, is only loosely constrained and the value of the [3 pa¬ 
rameter is unknown, a is known to strongly depend on the field 
strength, with a proportional to [3 1/2 (Hawley et al. 1995) or (3 1 
(Bodo et al. 2011) in the vertical field case. Similarly, we can ex¬ 
pect angular momentum transport to increase with the magnetic 
Reynolds number (Longaretti & Lesur 2010). The asymptotic 
values of the angular momentum rate we obtained are thus only 
valid for the f3 and Rm parameters we considered, and with a 
scaling that remains to be determined in the small Pm limit. 

In the case of the compressible simulations with a mean ver¬ 
tical field, we obtained a surprising ratio of Maxwell to Reynolds 
stress of about 1. Not withstanding any possible dependance 
with [3 and Rm, we noticed that this ratio is related to the com¬ 
pressibility of the fluid, as a more usual value of about 3 is ob¬ 
tained in the incompressible runs. It is also related to the bursty 
behavior of these two stresses: whereas a usual value is obtained 
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in the growing phase of the bursts, the Reynolds stress dominates 
during its decreasing parts, resulting in a mean R value of about 
1. This result is specific to the compressible simulations, but it 
implies that neither the Reynolds stress nor the Maxwell stress 
converge to the same value in compressible and incompressible 
simulations at low Pm. In fact, both stresses differ by about 50%, 
which we show can be attributed to their different behavior dur¬ 
ing the bursts. It is possible that the convergence of a at low Pm 
to the same value with this two types of flow is fortuitous. In 
order to solve that issue, higher Rm simulations are needed but 
are currently very costly. 

Next, we investigated the interest of using large eddy sim¬ 
ulations both in spectral and real spaces for such simulations. 
The simplest approach is to consider the implicit LES method. 
As expected, in all cases the smallest scales are not correctly 
handled, but it is possible to reproduce the large scales energy 
spectra by using a high enough resolution. To obtain an error 
smaller than 20% on a, the needed resolution corresponds to 
a fourth of the resolution needed when explicit viscosity is in¬ 
cluded. This corresponds to a decrease in CPU time by a factor 
256. To limit the computational cost of a simulation of a turbu¬ 
lent flow, explicit SGS models are also usually considered. The 
anisotropy of the rotating sheared flow and the turbulent cas¬ 
cade which differs from the Kolmogorov cascade, may indicate 
that simple SGS approaches are not well suited for MRI turbu¬ 
lence simulations. We tested the Chollet-Lesieur method which 
is used in spectral space and is applied directly on the spectrum 
components. This method is local in frequency space, the effec¬ 
tive viscosity being a function of the wave-number, and there is a 
limited inclusion of the backscatter. We found that good results 
are also obtained with this Chollet-Lesieur approach notwith¬ 
standing the anisotropy of the flow at small scales. 

Despite such positive results, a word of care is in order 
here. Although we showed that these methods are efficient to 
decrease the computational cost of MRI turbulence simula¬ 
tions, the needed resolution is still significant. For a magnetic 
Reynolds number Rm = 400 with a vertical mean magnetic field, 
the required resolution is 64 pts/H. This Rm is low compared to 
the values that are often chosen in the literature for which an 
even higher resolution is then necessary. To reach a turbulent 
flow dominated by Ohmic dissipation, a resolution of at least 
256 pts/H will be needed for Rm values of a few thousands. 
Moreover, we considered only two diagnostics, namely the an¬ 
gular momentum transport rate and the power spectra, to reach 
this conclusion and other diagnostics, such as helicity or corre¬ 
lation functions, were not considered. Such statistical quantities 
may well be incorrectly described in our LES for the resolu¬ 
tions we used. More generally, with such reduced resolutions, 
any small scale process is not correctly handled, and for instance 
the study of collision of grains in protoplanetary discs or mag¬ 
netic reconnection cannot be studied in such simulations. 

Overall we still conclude that LES can be used to limit the 
computational time of future simulations. Future works should 
account for density stratification, which will be particularly rel¬ 
evant in the context of global simulations. Large scale phenom¬ 
ena such as the ones identified in the ‘butterflies diagrams’ can 
strongly modify the flow properties and are susceptible to affect 
the LES approach. Dedicated simulations should be performed 
to quantify the interest of such methods in this case. 
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